# include "scene.h"
# include "discr.h"
# include <stdio.h>
# include <stdlib.h>

#ifdef _OMP_PRESS_BC_
# include <omp.h>
#endif

extern double inlet_velocity(double r, double t, short index);
extern double interpol_1d(FILE *lptr, double *x, double *u, int n, double xp, double *du);

# ifdef DENS_TDEP
extern double dens(double T);
# endif
# ifdef KVISC_TDEP
extern double kvisc(double T);
# endif
/*
# include <time.h>
# include <sys/timeb.h> */

extern void press_bc(const struct grid *g, short tflag, unsigned short iRK)
  {/* tflag = 0: current time step i.e. t; tflag = 1: next time-step i.e. t+dt. */
   short i, j, k, _k, nr, nz, nth, tid;
   double ***u, ***v, ***w, ***T, *r, *th, *z, _dwdt_iRK, L, density0, t, kvisc0, _z[PRESS_INTRP_N], (*_p)[PRESS_INTRP_N];
// struct timeb t1, t2;

// ftime(&t1);
#  ifdef TRACE_FILE
   fprintf(g->trptr, "Entering into 'press_bc()'...\n");
#  endif

   if( (iRK > 0 && tflag != 0) || (iRK == 0 && tflag != 1) || (iRK > 2) )
       {
        fprintf(g->lptr, "\nERROR: In 'press_bc()': (tflag, iRK) = (%d, %d).\nThis function is designed to handle only following cases:\n\t[tflag = 0, iRK = 1, 2] and [tflag = 1, iRK = 0]\n", tflag, iRK);
        exit(-1);
       }

   u = g->u;   v = g->v;   w = g->w; T = g->T; r = g->r;   th = g->th;   z = g->z;  nr = g->nr;   nz = g->nz, nth = g->nth;
   density0 = g->density0;
   kvisc0 = g->kvisc0;

   t = g->t + tflag*g->dt;

   _p = (double (*)[PRESS_INTRP_N])malloc(g->nCPU*PRESS_INTRP_N*sizeof(double));

   _z[0] = z[NzG];
   for(k = 1; k < PRESS_INTRP_N; k++) _z[k] = 0.5*(z[NzG-1+k]+z[NzG+k]);

   tid = 0;
#  ifdef _OMP_PRESS_BC_
   #pragma omp parallel private(i,j,k,_k,_dwdt_iRK,L,tid)
#  endif
     {
      k = NzG-1; /*** At the BOTTOM INTERFACE. It is used in evaluation of "L". ***/

   /* At the INLET: Numan BC */
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < g->nr_nozzle + NrG; i++) for(j = NthG; j < nth + NthG; j++)
           {
            switch(iRK)
               {
                case 0:
                   _dwdt_iRK = inlet_velocity(0.5*(r[i]+r[i+1]), t, 1);
                   break;

                case 1:
                   _dwdt_iRK = inlet_velocity(0.5*(r[i]+r[i+1]), t, 1)\
                               + C2*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 2);
                   break;

                case 2:
                   _dwdt_iRK = inlet_velocity(0.5*(r[i]+r[i+1]), t, 1)\
                               + (a31+a32)*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 2)\
                               + a32*C2*g->dt*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 3);
                   break;
               }

/*        if(iRK == 3) _dwdt_iRK = inlet_velocity(0.5*(r[i]+r[i+1]), t, 1)\
                             + (a41+a42+a43)*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 2)\
                             + (a42*C2+a43*a31+a43*a32)*g->dt*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 3)\
                             + a43*a32*C2*g->dt*g->dt*g->dt*inlet_velocity(0.5*(r[i]+r[i+1]), t, 4); */

            L =  0 - _dwdt_iRK - ZMOMT_TERM_3 + ZMOMT_VISC - ZMOMT_TERM_5;
            g->p[i][j][NzG-1] = g->p[i][j][NzG] - density0*(z[NzG+1]-z[NzG])*L;
           }


   /* At the SOLID-BASE: Numan BC */
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG+g->nr_nozzle; i < NrG+nr-g->nr_gap; i++) for(j = NthG; j < nth+NthG; j++)
            {
             L = ZMOMT_VISC - ZMOMT_TERM_5;
             g->p[i][j][NzG-1] = g->p[i][j][NzG] - density0*(z[NzG+1]-z[NzG])*L;
            }

   /* At the OUTLET: Dirichlet B.C. : It can be used because mass flow rate at the outlet is not specified. Because
      the boundary condition specified at the outlet is Numan BC and Dirichlet B.C. at the outlet can be specified. */
#     ifdef _OMP_PRESS_BC_
      tid = omp_get_thread_num();
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG+nr-g->nr_gap; i < NrG+nr; i++) for(j = NthG; j < nth+NthG; j++)
            {
             _p[tid][0] = P_ATM + density0*GRAVITY*h_WATER;
             for(_k = 1; _k < PRESS_INTRP_N; _k++) _p[tid][_k] = g->p[i][j][NzG-1+_k];
             g->p[i][j][NzG-1] = interpol_1d((FILE *)0, _z, _p[tid], PRESS_INTRP_N, 0.5*(z[NzG-1]+z[NzG]), (double *)0);
            }

   /* At the TOP-SURFACE: Numan BC. */
      k = NzG+nz-1; /*** Used in the evaluation of "L". ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided) nowait
#     endif
      for(i = NrG; i < nr+NrG; i++) for(j = NthG; j < nth+NthG; j++)
            {
             L = ZMOMT_VISC - ZMOMT_TERM_5;
             g->p[i][j][NzG+nz] = g->p[i][j][NzG+nz-1] + density0*(z[NzG+nz]-z[NzG+nz-1])*L;
            }


   /* At the SIDE-WALL: Numan BC. */
      i = NrG+nr-1; /*** Used in the evaluation of "L" ***/
#     ifdef _OMP_PRESS_BC_
      #pragma omp for schedule(guided)
#     endif
      for(k = NzG; k < NzG+nz; k++) for(j = NthG; j < nth+NthG; j++)
            {
             L = RMOMT_VISC;
             g->p[NrG+nr][j][k] = g->p[NrG+nr-1][j][k] + density0*(r[NrG+nr]-r[NrG+nr-1])*L;
            }
     }
   free(_p);

/* For ghost cells in theta-direction having possitive radius. */
#  ifdef _OMP_PRESS_BC_
   #pragma omp parallel for private(i,j,k) schedule(guided)
#  endif
   for(i = NrG; i < nr+2*NrG; i++) for(k = 0; k < nz+2*NzG; k++) for(j = 0; j < NthG; j++)
        { //Applying Symmetric boundary condition in theta-direction.
         g->p[i][NthG-j-1][k] = g->p[i][NthG+j][k];
         g->p[i][nth+NthG+j][k] = g->p[i][nth+NthG-1-j][k];
        }

/* For ghost cells having negative radius: applying symmetric B.C. across the center. */
#  ifdef _OMP_PRESS_BC_
   #pragma omp parallel for private(i,j,k) schedule(guided)
#  endif
   for(k = 0; k < nz + 2*NzG; k++) for(j = 0; j < 2*NthG + nth; j++) for(i = 0; i < NrG; i++)
          g->p[NrG-1-i][j][k] = g->p[NrG+i][j][k];

#  ifdef TRACE_FILE
   fprintf(g->trptr, "...'press_bc()' ends.\n");
#  endif

   return;
  }
